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Abstract 

The purpose of this paper is to examine the Lagrangian stochastic modehng of the fluid velocity 
seen by inertial particles in a non-homogeneous turbulent flow. A new Langevin-type model, com- 
patible with the transport equation of the drift velocity in the limits of low and high particle inertia, 
is derived. It is also shown that some previously proposed stochastic models are not compatible 
with this transport equation in the limit of high particle inertia. The drift and diffusion parameters 
of these stochastic differential equations are then estimated using DNS data. It is observed that, 
contrary to the conventional modeling, they are highly space-dependent and anisotropic. To inves- 
tigate the performance of the present stochastic model, a comparison is made with DNS data as well 
as with two different stochastic models. A good prediction of the first and second order statistical 
moments of the particle and fluid seen velocities is obtained with the three models considered. Even 
for some components of the triple particle velocity correlations, an acceptable accordance is noticed. 
The performance of the three different models mainly diverges for the particle concentration and 
the drift velocity. The proposed model is seen to be the only one which succeeds in predicting the 
good evolution of these latter statistical quantities for the range of particle inertia studied. 
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I. INTRODUCTION 



In the past years, several numerical methods have been employed to study the dispersion 
of solid particles in turbulent flows. Generally, small enough particles are considered in order 
to treat them as point-particles.-'^'^'^ Assuming that the drag force is only of importance, the 
link between the motion of an inertial particle and the carrier fluid is given by the following 
system of equations : 



(1) 
(2) 



dt 

dvp.i Ui - fp. 



dt Tp 

where Xp j and the particle position and velocity, Tp is the particle relaxation time 

which is expressed in terms of the drag coefficient and of the magnitude of the relative 
velocity, and Ui = Ui{'Xp,t) is the fluid velocity at the particle location. Under these consid- 
erations, the main difficulty then lies in the proper computation of the fluid velocity at each 
particle location. The first possibility is to use Direct Numerical Simulation (DNS).-"^ This 
technique gives the best estimation of the fiuid velocity seen by particles. Nevertheless, it ne- 
cessitates very high computational ressources. A more affordable numerical way is provided 
by Large Eddy Simulation (LES).- Contrary to DNS, a model which takes the residual fluid 
dynamic (i.e. at the sub-grid scale) into account should be used to predict the instantaneous 
fluid velocity seen by particles. Finally, when the computational cost of this latter technique 
is still too high, it is then possible to make use of macroscopic numerical simulation such 
as Reynolds- Averaged Navier-Stokes (RANS). The use of a RANS-Lagrangian method to 
describe the motion of solid particles in a turbulent two-phase flow necessitates to generate 
the fluctuating velocity of the carrier phase at particle location.- In this framework, aver- 
aged quantities such as the mean velocity and some of the mean turbulent characteristics of 
the carrier phase are determined by solving the Reynolds- Averaged Navier-Stokes equations. 
The time integration of the equations governing the motion of inertial particles [Eqs. ([T]) and 
([2])] requires the knowledge of the instantaneous velocity of the fluid at the particle location. 
The reconstruction of the random nature of the fluctuations along inertial particle trajecto- 
ries can be achieved, for instance, using a stochastic Lagrangian models.-'^'^'^'^'^ Most of 
these models for the simulation of turbulent two-phase flows involves speciflc formulations 
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based on the Langevin model which can be written in a general form as 



dui = Aidt + BijdWj , (3) 

for the instantaneous fluid velocity at particle location. In this latter stochastic differential 
equation (SDE) , Ai is the drift vector, Bij is the diffusion matrix, and dWj are the increments 
of a vector- valued Wiener process with independent components. Some important properties 
of these increments are that they are non-differentiable and normally distributed with mean 
{Wi{t + dt) - Wj{t)) = and variance {[Wi{t + dt) - Wj{t)f) = dt6ij.-^ In order to predict 
the fluid velocity seen by inertial particles, one has to model the drift vector and the diff'usion 
matrix. In the present study, we focus on the models for the drift vector proposed by Simonin 
et alA 

and by Minier and Peirano^ 

In these expressions, u is the kinematic viscosity, pf is the fluid density, p stands for the 
pressure, Ui is the fluid velocity, Vp^i is the solid particle velocity and Gij is the drift matrix. 
The difference between both models lies in the form of the third term which describes the 
crossing-trajectory effect.— In the model proposed by Simonin et al.—, this term is written 
as a function of the instantaneous relative velocity between the particle and the fluid while 
Minier and Peirano^ suggested to express it as a function of the mean relative velocity. 
It has to be noted that in the limit of low particle inertia (r^ ^ 1), both models give the well- 
known Generalized Langevin Model (GLM) derived by Pope^ to predict the motion of fluid 
particle in a turbulent flow. Nevertheless, as it will be shown in Sec. II, these two previous 
stochastic models are not compatible with the transport equation of the drift velocity (mean 
fluctuating fluid velocity at particle location) for large particle inertia. In order to correct 
this discrepancy, a new form of the drift vector is proposed. In Sec. Ill, the method used 
to derive the drift and difi'usion parameters of these stochastic models is described, and the 
estimated values obtained using DNS data are presented. The performance of the proposed 
functional form of the drift vector is then assessed by comparison with DNS data in Sec. 
IV. Finally, concluding remarks are given in the last section. 
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II. EXACT AND MODELED TRANSPORT EQUATIONS OF THE DRIFT VE- 
LOCITY 

In this section, we study the Langevin models proposed by Simonin et al— and Minier 
and Peirano^ through the transport equation of the drift velocity in the limit of low and 
high particle inertia. Based on this study, a new model for the drift vector Ai is proposed. 



A. Degenerate equations for low and high particle inertia 

Let us consider the gas-solid flow from an Eulerian (macroscopic) point of view. The 
exact transport equation for the statistical moments of the particle and fluid seen velocities, 
as well as for the fluid seen-particle velocity correlations, can be derived, for example, from 
the transport equation of a joint probability density function for the particle and fluid seen 
velocities.— '^'^'^ The exact transport equation of the drift velocity, (£t^) = (£tj — (uj)), can 
be written as^ 



'Dt \_dx, 



d_ 



Wp,j) — i^ppp) 



- appp {{vpj) - {uj)) ^— - appp {uj) 

f I d{p) 92 {ui) \ /duA 

where Dp{-)/ Dt = d{-)/dt + {vpj) d{-)/dxj, and ap is the particle volume fraction. 
In the limit of vanishing particle inertia (rp <^ 1), a solid particle behaves like a fluid particle 
tracer. Its velocity is equal to that of a fluid particle, the statistical moments of the fluid 
and particle velocities are thus identical. Moreover, the drift velocity is zero since this kind 
of particles samples homogeneously the turbulent flow fleld and the particle volume fraction 
is constant if the particles are uniformly distributed initially. As a consequence, it can be 
found, from equation ([6]), that the average of the time variation of the fluid velocity seen 
becomes equal to 

dui \ _ Id (p) ^ d"^ {uj) 



dt I pf dxi dxjdxj 
The averaged Navier-Stokes equations are thus recovered. 

The opposite limit case, i.e. high particle inertia (r^ ^ 1), is also of great importance when 
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studying gas-solid flows since the trajectories of such particles become completely indepen- 
dent of the fluid motion. In such a case, the particle velocity remains nearly identical to its 
initial value, the fluid seen-particle velocity correlations, the second and higher statistical 
moments of the particle velocity as well as the drift velocity tend to zero. Moreover, the 
particle volume fraction keeps a constant value across the fluid flow if particles are initially 
distributed uniformly. In the present study and without loss of generality, the velocity of 
these high inertia particles is considered identical to the mean fluid velocity. Under these 
considerations, it can be found from equation ((6l) that 



These two previous equations give the asymptotic limits of the average time derivative of 
the fluid velocity seen by particles. Besides, they can be used in order to verify that the 
Langevin models generally used to predict the evolution in time of this fluid velocity are 
correct in the limits of low and high particle inertia. 

For example, let us consider first the model of the drift vector proposed by Simonin et alA, 
i.e. equation dH). The average of equation ([3]) in the limit of low particle inertia yields 



since (BijdWj) = 0. Thus, the model of the drift vector by Simonin et al— is able to produce 
the correct limit of the mean time increment of the fluid velocity seen in this particular case. 
When the particle inertia becomes high, the same expression is obtained from this model. 
In comparison with the exact one given by equation ([8]), it is noticed that the divergence 
of the Reynolds stress tensor is missing. Therefore, we can expect that some discrepancies 
could occur in the prediction of the momentum exchange between the dispersed and carrier 
phases for high particle inertia. Considering the model proposed by Minier and Peirano^, 
it can be seen that the expressions of the mean time increment of the fluid velocity seen in 
the limits of low and high particle inertia are identical to those obtained with the model 
of Simonin et al.— This model is thus also not compatible with the transport equation of 
the drift velocity in the limit of high particle inertia. In the next section, a new model of 
the drift vector which makes possible the prediction of the theoretical limits given above is 
proposed. 






1 d{p) 
Pf dxi 



+ y 



dxjdxj ' 



(9) 
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B. Proposal of a new model 



From the conclusions drawn in the previous section, we have designed a new model which 
gives the proper limits of the drift velocity transport equations in the cases of low and high 
particle inertia. This model for the drift vector Ai is 

(10) 

This drift vector is mainly different from those proposed by Simonin et ai— and Minier and 
Peirano^ due to the presence of the term d {{u'^'p^i^) — i'^'i'^'k)) I^Xk (i.e. the divergence of 
the difference between the fluid-particle covariances and the Reynolds stresses). 

First, it has to be noted that in the limit of low particle inertia the model proposed is 
also identical to the GLM derived by Pope^ because (u'^Vp^f^) (u[u'j}i. In addition, the 
modeled averaged time increment of the fluid velocity seen using Eq. (fTO]) has the proper 
limits since 

/ dui \ _ 1 d{p) ^ d'^ {ui) 

when Tp <C 1, and 



\ dt ) pf dxi ~^ dxjdxj ' ^ ^ 



dt I Pf dxi dxjdxj dxj ^ * ^' ' 



for Tp > 1. 



The introduction of a supplementary term, which is a function of the Reynolds stresses, in the 
drift vector has been motivated by the necessity for the stochastic model to be consistent with 
the transport equation of the drift velocity in the limit of high particle inertia. Nevertheless, 
this additional term had to vanish in limit of low particle inertia in order to keep the model 
similar to the GLM. This has naturally led us to add the fluid-particle covariance tensor 
which tends to the Reynolds stresses when Tp ^ 1 and to zero when Tp ^ 1. Moreover, 
incorporating the proposed model in the transport equation of the drift equation, it can be 
seen that this new term moment is physically consistent with the others. 

At this point, we would like to emphasize the fact that the present model should be 
more suitable for predicting the fluid velocity seen by large solid particles than the models 
proposed by Simonin et al— and Minier and Peirano^, however, the presence of the fluid- 
particle covariances in the expression increases the degree of complexity of the stochastic 
model. 
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It is also worth mentioning that the SDE for the fluctuating fluid velocity seen derived 
from this model (this SDE is presented hereafter) presents similarities with the one proposed 
recently by Bocksell and Loth—. In this latter study, they concluded that a "drift correction", 
which is a function of the fluid seen and particle velocities, should be included in the SDE in 
order to correctly predict the concentration proflles of flnite-inertia particles in a turbulent 
boundary layer. In fact, the last term in Eq. (|TOl) plays this role. 

Before evaluating the performance of the proposed model, the procedure used to specify 
the parameters of the stochastic equation, i.e. the drift and diffusion matrices {Gij and Bij), 
is presented. 

III. DETERMINATION OF THE DRIFT AND DIFFUSION MATRICES 
A. Theoretical formalism and assumptions 

In order to test the capability of the proposed form for the drift vector to model the 
turbulence seen by inertial particles, the values of the components of the drift and diffusion 
matrices, Gij and Bij, have to be specifled. In stationary homogeneous isotropic turbulence 
and without a mean relative motion between the dispersed and carrier phases, the drift term 
is modeled has the inverse of the integral time scale of the fluid seen. The diffusion matrix 
is generally supposed independent of the particle inertia and is expressed as a function of 
the Kolmogorov's constant and dissipation rate of the mean turbulent kinetic energy accord- 
ing to the Kolmogorov similarity theory for the second-order Lagrangian velocity structure 
function in the inertial subrange.— It has to be noted that this model for the diffusion term 
is strictly valid in the limit of vanishing particle inertia and for high Reynolds number tur- 
bulent fluid flows. 

In the case of non-homogeneous turbulence, the speciflcation of the drift and diffusion matri- 
ces is even more complex. There are no models for these quantities which take properly the 
properties of such a turbulence into account. Consequently, we propose in the present study 
to determine Gij and Bij using data extracted from our channel flow DNS computation. 
A similar method to the one proposed in the study by Pope^, which was devoted to the 
prediction of fluid particle trajectories in a turbulent homogeneous shear flow, is followed. 
Firstly, in order to apply this method, the stochastic differential equation for the fluctu- 
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ating fluid velocity at the solid particle location has to be derived from Eq. ([3]). Since 

du[ = dui — d (ui), it can be shown that 

du\ = Giju'jdt + BijdWj + ^' dt , (13) 



when the model by Simonin et al— is used while 



du[ = G.jUjdt + BijdWj + V dt , (14) 

UXj 

with the present model for Ai [Eq. (fTOl)]. In these equations, Gij = Gij — d {u-i) jdxj. The 
model of Minier and Peirano^ will not be considered in the rest of the present study for two 
reasons. This model suffers from the same drawback in the limit of high particle inertia as the 
one suggested by Simonin et al.—, consequently, only one of these models can be examined. 
In addition, their model was designed to be used for the prediction of the instantaneous fluid 
velocity seen by particles and is thus not of practical use for predicting the fluctuating part. 
Secondly, since the method proposed by Pope^ is strictly valid for homogeneous turbulent 
flows, an assumption has to be made in our case. We will assume that the turbulence is 
locally homogeneous so that the spatial derivatives of the turbulent statistics vanish. This is 
a strong assumption, however, we are interested in a fairly good and simple approximation of 
the drift and diffusion matrices in order to test stochastic models. As far as we know, there 
is no other simple method to determine the parameters of this particular type of stochastic 
models due to the non-homogeneity of the turbulent flow studied. Moreover, it will be 
shown later from the stochastic simulations of the gas-solid flow that this approximation 
leads to very good results. Besides, it should be also noted that it is under this assumption 
that Walpot et al— recently derived the drift matrix of a stochastic equation predicting 
the fluctuating velocity of fluid particles for a turbulent pipe flow. Nonetheless, it has to be 
mentioned that other stochastic models, motivated by the works of Wilson et al.—, Durbin— , 
and Thomson—, have been suggested to tackle the problem induced by the non-homogeneity. 
More details can be found in Iliopoulos and Hanratty— , Iliopoulos et al.—, and references 
within. 

Assuming the turbulence as locally homogeneous, the drift matrix can be expressed from 
Eq. ^ or Eq. ^ as 

G = -{7T\ (15) 
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where (■)^ denotes the transpose and Tij is the matrix of the decorrelation time scales of 
the fluid seen which is deflned as 

poo 

'^^J= / {U'^',)-' (Kmrit)) dt , (16) 

with {u[u'i.)~^ being the i — k component of the inverse of (uu^). In order to obtain the 
drift matrix, 7ij has been computed from DNS data. 

To determine the diffusion matrix, we have to consider the transport equation of the sec- 
ond order statistical moment of the fluid velocity seen by particles which is described by 
equation ([3]). This transport equation, which can be derived from the transport equation 
of the joint probability density function for the particle and fluid seen velocities,— has the 
following form 

^ ^ ' y dxk pf dxj dxkdxk 
'P 1' \ dxk Pf dxi dxkdxkj 



dxk ^' dxk 

+ apPp {Aiu'j + Aju[) + apPp {BikBjk) . (17) 

Note that this transport equation is generally written for convenience in terms of a fluctu- 
ating fluid velocity seen deflned as u'- = Ui— {ui) while we deflned it as u'^ = Ui— {u.j) in the 
present study. The other form of the transport equation can be thus found by introducing 
the relation u'^ = u'l + (n^) in Eq. (flTll . 

Let us now write the drift vector, Aj, in a compact form as 

When Cj = 0, the model proposed by Simonin et al— is recovered while the new proposed 
model is obtained if Cj = 9 {{u[v'^ f}i — (m^m'^)) /dxk- Introducing the expression of the drift 
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vector in Eq. ffTTll yields 




dxk 



(19) 



The gas-solid channel flow being statistically stationary and homogeneous in the streamwise 
and spanwise directions, the difl'usion matrix can be expressed, under the local homogeneity 
assumption, from Eq. (fT9l ) as a function of the drift matrix 



Here, it should be noted that Bfj does not determine uniquely Bij. Nevertheless, Bfj 
will produce a unique set of statistical moments of the fluctuating fluid velocity seen by 



The parameters of the Langevin model can be thus expressed as a function of the decor- 
relation time scales and second order statistical moment of the fluid seen by particles. 

B. Results 

In order to evaluate the parameters of the Langevin model, data extracted from a direct 
numerical simulation of a gas-solid channel flow have been used. The direct numerical 
simulation was conducted at a Reynolds number Rcb = 2280 (based on channel half-height S 
and bulk velocity Ub) corresponding to a Reynolds number based on the wall-shear velocity 
(ur) equals to Rcr ~ 155. The results used in the present study are coming from the same 
numerical computations presented in Marchioli et al.— to test the prediction of particle 
dispersion by different DNS codes. Therefore, only the main characteristics of the gas-solid 
flow simulation are given here. The numerical simulation of solid particle trajectories was 
restricted to spherical particles smaller than the smallest turbulent length scale. Therefore, 
we made use of the point-force approximation. In the present study, the particle-particle 
interactions as well as the turbulence modulation were disregarded (one-way coupling). In 
addition, the added mass, history and lift forces were neglected in the particle equation 
of motion since the ratio between the particle and fluid density obeys Pp/ Pf ^1. In the 




(20) 



particles 



Therefore, we suppose in this study that Bij is symmetric. 
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present study, only the non-linear drag force, estimated from the correlation of Morsi and 
Alexander—, was considered. 

Simulations were run for three sets of particles characterized by different Stokes particle 
response times in wall units, r+ = 1, 5 and 25 [quantities in wall units are normalized 
with the viscous scales (i.e. the wall-shear velocity Ur and the viscous lengthscale z//u^) 
and indicated by the superscript (■)^]. The corresponding dimensionless diameters were 
dp/6 = 1 X 10~^, 2.2 X 10~^ and 5 x 10~^, and the density ratio was equal to Pp/pf = 1000/1.3 
for the three sets. Statistics on the dispersed phase were started after a time lag necessary 
for particle statistics to reach a stationary state. 

In figures [H and O the components of the drift matrix in wall units are plotted as 
a function of the wall-normal coordinate and for the three different particle inertia. 
Before commenting the results, we have to emphasize that not too much attention should 
be paid to the behavior of the Langevin model parameters near the wall since they were 
derived under the assumption of local homogeneity. This approximation is certainly not 
correct in this region of strong gradients of the turbulent statistical moments. From the 
diagonal components of Gf^ shown in Fig. [H it can be observed that the magnitude of G22 
and Gj^3 decreases monotically with increasing y~^. Moreover, the particle inertia is shown 
to not have a significant effect on these components. These trends are quite different for 
Gfi since its magnitude is seen to have a local minimum at y''" ~ 10 whatever the particle 
inertia, and then for y"*" > 40, it decreases with increasing Concerning the non-diagonal 
components of Gfj plotted in Fig. [21 we note that Gf2 is maximum in the near- wall region 
and tends to zero at the channel center. The particle inertia has a quite important effect 
on that component while G21 is zero across the channel whatever the particle inertia. 
Finally, it has to be mentioned that our estimation of the drift matrix for the lowest particle 
inertia is qualitatively in good agreement with the results obtained by Walpot et al— in 
their study of a Langevin model for predicting the fiuctuating velocity of fiuid particles in 
a turbulent pipe fiow. 

The results obtained for the diffusion matrix are given in the form of Bf- = BikBjk in 
figure [3l Contrary to the conventional modeling assumption, it is found that the diffusion 
term Bf, is significantly anisotropic for < 100. A similar observation was previously made 
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by Pop©^ in a study of the stochastic Lagrangian modeling of fluid particle trajectories in 
a homogeneous turbulent shear flow. Besides, the results show that the components of Bf- 
tend towards zero close to the wall and have a maximum located approximately at y"*" = 25. 
Consequently, according to these results, Bf- cannot be modeled as a function of kinetic 
energy dissipation rate moderated by a constant since the dissipation of the kinetic energy 
is maximum at the wall. Concerning the inertia effect, it is observed that the values of Bl^ 
are identical for 7^=1 and 5 and increase when r+ = 25. This is not the case for the 
components iJfg and i^fg since similar results are obtained for r^*" = 5 and 25 particles while 
the magnitude of these components is higher for the lowest particle inertia. Regarding the 
non-diagonal component, the particle inertia effect is seen to only change its minimum. Using 
these results for the drift and diffusion matrices, the performance of the present stochastic 
model is examined in the next section. 

IV. EVALUATION OF THE STOCHASTIC MODELS 
A. Presentation of the test 

To assess the performance of the proposed Langevin model, a comparison between re- 
sults obtained from a stochastic simulation and those extracted from the direct numerical 
simulation of a gas-solid channel flow has been conducted. The stochastic simulations have 
been carried out for three different forms of the drift vector in order to investigate the effects 
of this term on the predicted dispersed phase statistics. The expressions of the stochastic 
differential equation corresponding to these models can be put in the following compact 
form: 

du'i = Giju'jdt + BijdWj + Didt . (21) 

Consequently, Di = d {u[u'i^ / dxk gives the SDE obtained using the model proposed by 
Simonin et alA [Eq. (fT3l) ]. = d {u'^v'p^i^') /dxk is the second form which is derived using the 
proposed model for the drift vector [Eq. (fT4ll ]. and = corresponds to the third model 
considered. This last form is less cumbersome than the two others and does not need to 
know beforehand the fluid Reynolds stresses or the fluid-particle covariances. In addition, 
it should be noted that this model only gives the proper limit of the average of the time 
variation of the fluid velocity seen when Tp ^ 1 [see Eq. ([HD]. 
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B. Stochastic simulation 

For the stochastic simulation of the gas-sohd flow, the mean fluid motion was calculated 
by means of a Reynolds- Averaged Navier-Stokes model. Closure of the Reynolds stresses 
is achieved using the Non Linear Eddy Viscosity Model,— so that turbulence anisotropy 
is taken into account. In order to have a better precision in the near- wall region where 
the viscosity effects have to be taken into account, modifications of the standard k — e 
model following the recommendations of Myong and Kasagi— have been made. This model 
introduces "damping" functions that allow the transport equations of k and e to be valid 
close to the wall. A complete description of the RANS model used here is given in the work 
of Carlier et al.—. The mean fiuid velocity, (lii)^, predicted by this model is compared in 
figure [4]to the one obtained by DNS. Despite a slight overprediction of {ui)'^ by the RANS 
model in the logarithmic region, the results can be considered to be in good accordance. 

After the Eulerian computation of the mean fiuid velocity, the stochastic Lagrangian 
simulations have been performed by tracking as much as 2.10^ solid particles in order to get 
enough statistical information in each cell of the domain to calculate the mean dispersed 
phase statistics. The particle characteristics as well as the equation of motion used for 
calculating the trajectories are identical to those of the DNS computation. 
Nevertheless, contrary to the gas-solid DNS which has been conducted in a bi-periodic 
domain, a finite streamwise length channel has been considered for the stochastic simulation. 
This length has been chosen to be equal to 10 m (~ 500 times the channel half-width) in 
order to obtain dispersed statistics, calculated at the outlet, which are independent of the 
distance to the inlet. 

The fiuid velocity fiuctuation at the particle location has been determined by integrating in 
time the stochastic equation [Eq. (|2T]1 ] in a semi-analytical way. The stochastic part is firstly 
disregarded and the time increment of the velocity can be thus analytically obtained from 
the resulting system of coupled equations. The stochastic term increment is then estimated 
using an Euler scheme and added to the analytical solution. The simulations have been 
carried out using a time step being equal to rp/25 for r + = 5 and 25. For the r + = 1 
particles, the time step was chosen to be rp/5 in order to limit the computational cost. One 
should note that during these stochastic simulations, the time step is always lower or of the 
order of the smallest velocity timescale characterizing the present fiow. This choice is also in 
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accordance with the guideline given by Sommerfeld— . The values of Gij and Bij have been 
linearly interpolated at the solid particle location from the data presented in the previous 
section. These coefficients being unknown at the wall, a linear extrapolation has been chosen 
to estimate them near the walls (0 < i/+ < 3.1). In addition, the mean turbulent statistics, 
which appear in the three tested stochastic models, have been also calculated at the particle 
position using a linear interpolation of data extracted from the DNS computation in order 
to not introduce additional modeling uncertainties. 

C. Numerical results 

The first result we present is the concentration (more precisely the number density) of 
solid particles across the channel width. The DNS and stochastic simulations conducted 
with the three different models are compared in figures [5](a-c) for r+ = 1, 5 and 25 
respectively. These results are interesting since they reveal that the drift model has an 
important effect on the particle distribution in the channel. The DNS data show that 
the particle concentration increases with increasing inertia (in the particle inertia range 
studied). This behavior can be seen as a preferential concentration effect at the macroscopic 
scale. Of course, it is different from the the local effect which occurs at smaller length 
scales .-'^'^'^ 

The comparison with the stochastic simulation shows also that the model of the drift 
vector proposed by Simonin et alA [Eq. (l2Tl) with Di = d {u[u'j^) /dxk] is not able to 
reproduce the accumulation of the larger particles in the low turbulent intensity regions. 
The concentration profiles remain quite uniform whatever the particle inertia. It seems 
that the presence in the SDE of the term Di = d {u[u'j^) /dxk (which does not vary as 
a function of particle inertia) prevents the accumulation. On contrary, the model with 
-Dj = produces a segregation of the particles in the near-wall region even for the lowest 
particle inertia. This is in contradiction with the law of conservation of mass since these 
particles, which can be assimilated to fiuid particle tracers, have to approximately be 
uniformly distributed.— This non-physical behavior, which is called spurious drift effect, 
was observed by Wilson et al— and Maclnnes and Bracco^ in stochastic Lagrangian 
simulations of tracer particles in non-homogeneous turbulent fiow. In fact, it is due to an 
inconsistency between the stochastic Lagrangian model and the Navier-Stokes equations 
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which causes a misrepresentation of the averaged time derivative of the fluctuating fluid 
particle velocity. A more detailed presentation of this effect can be found in Popo^, 
Thomson—, Guingo and Minier— , and references within. Despite this major problem for 
low particle inertia, this stochastic model predicts reasonably well the concentration of 
the r+ = 25 particles. This confirms that this latter model should be more appropriate 
to estimate the fluid velocity seen by large particle inertia. Finally, it is noted that the 
results obtained with the present model for r + = 1 and 25 are in good accordance with 
the DNS data while important discrepancies are observed for r+ = 5 when < 2. 
Nevertheless, the main point is that this model reproduces qualitatively quite well the effect 
of inertia on the particle concentration. There is no spurious drift effect for low inertia 
and an increase of the concentration in the near- wall region is noted for the higher particle 
inertia. This result is of importance since the correct prediction of particle flux in wall- 
bounded turbulent flows is decisive when studying the complex process of particle deposition. 

The first order statistical moment of the particle velocity is plotted in figures [6](a-c). 
From Fig. [6](a), it is seen that the mean velocity of the smallest particle inertia studied 
is quite well predicted across the channel by the different models of the drift vector. 
Nevertheless, this velocity is overestimated in the viscous sublayer for < 2. The observed 
discrepancy can be reasonably attributed to the combination of two approximations. The 
first one concerns the local homogeneity assumption made to derive the parameters of the 
stochastic which does not hold in this region. The second one is the linear extrapolation 
used to estimate these parameters at the particle location. These remarks should be kept in 
mind throughout the presentation of the results. In the buffer and logarithmic regions, the 
present model as well as the one of Simonin et alA give similar results which are slightly 
greater than those of the DNS. This is due to the fact that the RANS model slightly 
overestimates the mean fluid velocity given by the DNS. Concerning the results obtained 
with the third model tested, i.e. Eq. (12X1) with Di = 0, we note that they are similar 
to those of the two other models except in the buffer region where the mean particle is 
lower. It is believed that this interesting difference is due to the fact that this model is 
not compatible with the transport equation of the drift velocity in the limit of low particle 
inertia. This incompatibility certainly gives rise to a wrong estimation of the drift velocity 
which should be quite low for this kind of particles. Since the mean particle velocity is 
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lower than expected, it can concluded that this model generates fluctuations of the fluid 
velocity whose average is negative. 

The mean velocity of the r + = 5 particles is shown in Fig. [6](b). Firstly, the present model 
and the one with Di = give identical results. The mean particle velocity is also slightly 
overestimated in the buffer and logarithmic regions. Nonetheless, the discrepancy observed 
for the smaller particles in the near-wall region is attenuated. This is certainly due to the 
less important sensitivity of the r + = 5 particles to the fluctuating fluid velocity. Secondly, 
it is noticed that the model of Simonin et alA causes a too high particle velocity in the 
buffer region. The incompatibility of this model with the transport equation of the drift 
velocity in the limit of large particle inertia can be invoked. Nevertheless, this explanation 
has to be taken with caution since one could argue that Tp = 5 particles do not belong to 
the category of large particles. At this stage of the study, no clear conclusion can be drawn. 
The mean velocity of r+ = 25 particles is plotted in Fig. [6](c), the three models give 
approximately the same results in the buffer and logarithmic regions. The predicted 
velocity is once more slightly higher than the one extracted from the DNS. In the viscous 
sublayer, this overprediction is also observed for the model of Simonin et al— while the 
results with the two other models are in a good accordance with the DNS. 

In order to better understand the influence of this three different stochastic models, the 
particle velocity root mean square (rms) has also been computed. The results obtained are 
shown in figures WM The estimated values of the streamwise particle velocity rms for the 
Tp = 1 particles [Fig. [7](a)] are seen to be in very good accordance with the DNS results. 
For the r + = 5 particles [Fig. [7|(b)], some differences are noted for < 30. Surprisingly, 
better results are obtained with the simplest model studied [i.e. Eq. fl2Tl) with Di = 0]. 
The two other models overestimate the streamwise particle velocity rms. The difference 
is roughly of the order of 5%. The results obtained for Tp = 25 [Fig. [7|(c)] show that the 
present model and the one with Di = predict accurately this particle velocity statistical 
moment. Some important discrepancies are noticed with the model by Simonin et al— for 

< 15. Concerning the wall-normal component plotted in figure [HI it can be seen that 
the accordance with the DNS results is very good whatever the model used. The predicted 
values are slightly lower than those given by the DNS for = 1 and 5, however, this 
difference is not significant. The inertia filtering effect is consequently well reproduced by 
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the stochastic models. For sake of conciseness and due to its strong qualitative similarity 
with the wall-normal component, the spanwise velocity rms is not shown. The non-diagonal 
component of the particle velocity rms is presented in Fig. [9l More important differences 
are observed for this component than for the others. Concerning the smallest particle 
inertia, the three models are in good agreement with the DNS data. For higher inertia, 
the results obtained from these models diverge for 10 < y'^ < 50. In this region, the 
model by Simonin et alA underestimates the magnitude of the minimum of particle kinetic 
shear stress {vp^iv'^^^. This trend is also obtained with the two other models, however, the 
difference is less. It should be also noted that the magnitude of (^^p ifp 2) given by DNS 
is higher than that of the three stochastic models for y+ > 50. A possible reason for this 
disagreement is that the wall-shear velocity used to normalized the quantities shown is 
not perfectly identical in the DNS and stochastic simulations. Finally, it can be remarked 
that the differences noticed previously are less for the = 25 particle except for y"^ > 50 
where the underestimation of the magnitude particle kinetic shear stress is of the same order. 

To complete this comparison of the particle velocity statistics, the prediction of a 
higher statical moment (the triple particle velocity correlations {v'^^jV^^v'^ j^) has been also 
investigated. This will give an idea of the capability of stochastic modeling. The agreement 
with the DNS data is not expected to be as good as for the other statistical moments shown 
due to the assumptions made for the derivation of the parameters of the stochastic models. 
The streamwise triple particle velocity correlation is plotted in figure [lOl Surprisingly, 
the three stochastic models are seen to be able to predict very well this correlation for 
y'^ > 30 and whatever the particle inertia. Moreover, the models well estimate the inertia 
effect since the obtained maximum of this correlation (located at ~ 10) increases with 
increasing inertia. Nonetheless, the magnitude of this maximum is not well predicted since 
the present model as well as the one by Simonin et al- clearly overestimate it while an 
underestimation is noted for the model with Di = 0. From figure \T2i it can be observed 
that {v'p 2v'p iVp i') is also quite well predicted by the present model and the one with Di = 
across the channel and whatever the particle inertia. Concerning the model of Simonin 
et alA, more important discrepancies arise for y~^ < 20 and = 5 and 25. For the 
correlations (^p,2^p,2^p,2) (^'p,2^p,i^'p,2)i shown in figures [TT] and [T3l a similar trend is 
noted. First, the results given by the three stochastic models are almost identical. Secondly, 



17 



the agreement with the DNS data is very good when < 20 whereas the magnitude 
of these two correlations is significantly underestimated in the rest of the channel. The 
last non-zero correlation computed is (^p,2^p,3^p,3) (Fig- El)- There are major differences 
between the DNS and stochastic simulations. The predicted correlation is generally smaller 
than the one extracted from the DNS. In addition, for < 40 and the smallest particle 
inertia, the wrong sign of (^p.2^p,3^p.3) is given by the stochastic models. Nonetheless, the 
results obtained with these models are in a good qualitative agreement for y^ > 40, and 
the inertia effect, which causes a decrease of this correlation, is well estimated. 

The second part of this comparison between the DNS and stochastic simulations is 
devoted to the statistics of the fluid seen velocity. The first statistical moment studied is 
the drift velocity (n^). In figures fTSl and [T6l the streamwise and wall-normal components 
are presented. A better prediction of the drift velocity from the present model is expected 
since it has been previously shown that it is the only model of the three considered which is 
compatible with the transport equation of this velocity in the limits of low and high particle 
inertia. From the results obtained for the streamwise component, it can be seen that the 
three models are unable to correctly estimate it for the = 1 particles when y^ < 40. 
The present model and the one of Simonin et alA predict quite correctly {u[) in the rest of 
the channel while significant difference with the DNS data are still noted using the model 
with Di = 0. This is in line with the fact that this latter model is not compatible with the 
transport equation of the drift velocity when ^ 0. As shown in the first part of this 
study, the two other models become identical in this limit. Nevertheless, there are some 
differences. This is due to the value of the particle inertia studied which is not enough low to 
observe the convergence of these two models. It is confirmed by the DNS data since the drift 
velocity of the r + = 1 particles is non-null while this velocity has to vanish when r+ — 0. 
From figure [TSlfb). it is apparent that the model with Di = better predicts the drift 
velocity for this kind of particles, however, there are important inconsistencies in a large 
part of the channel. The results given by the present model are in satisfactory agreement 
with the DNS data, whereas important qualitative and quantitative discrepancies are noted 
for the model by Simonin et al- when 10 < y"*" < 30. In the case of the highest particle 
inertia, this latter model and the present one give surprisingly similar results which are in 
acceptable agreement with the DNS data while the predictions by the simple model with 
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Di = are poor. It could have been expected that this latter model would lead to better 
results for these particles. Nonetheless, the inertia of the particle studied is not enough high 
to show that the present model and the one with Di = should predict the drift velocity 
more accurately than the model by Simonin et alA. The capability of the three models can 
be better distinguished from the results of the wall-normal component of the drift velocity. 
The stochastic model with Di = predicts a null drift velocity whatever the distance to the 
wall and the particle inertia. It is in complete disagreement with the DNS data. The results 
obtained with the model by Simonin et al.— are in satisfactory accordance. Nevertheless, 
important discrepancies begin to arise as the particle inertia increases. Contrary to these 
two models, the expression proposed to estimate the drift vector of the stochastic equation 
leads to a good estimation of this drift velocity whatever the particle inertia. These 
observations can explain the predictions of the particle concentration by the three models 
making use of simple physical considerations. The wall-normal drift velocity given by the 
stochastic model with = is null. It has been also shown that this model predicts an 
increase of particle concentration in the near-wall region whatever the particle inertia. This 
non-physical increase in the case of the lowest particle inertia is induced by the fact that 
there is no mean force to counteract the accumulation of the particles in low-turbulence 
regions (the so-called turbophoresis effect^i^). To explain the uniform concentration 
obtained with the model by Simonin et al.— whereas an increase should be observed near 
the wall in the case of the = b and 25 particles, similar considerations can be put 
forward. The predicted drift velocity given by this model, which is higher than that of 
the present model, is certainly too high to make possible the accumulation of these particles. 

The last statistical moment considered in this evaluation of the proposed stochastic model 
is the fluid seen-particle velocity correlations (^i^pj)- Concerning the diagonal components 
presented in Figs. [T7l and [THl the agreement of the stochastic simulations with the DNS data 
exhibits the same trends as for the particle velocity rms. The three models almost perfectly 
reproduce the inertia effect on the wall-normal and spanwise components. It should be 
noted that due to the strong similarity between these two components, only the results 
obtained for {u^v'^^^) ^^e shown (see Fig. [181). From figure [TTl it is noted that the streamwise 
component is in good accordance with the DNS data in the case of the smaller particle 
inertia. Differences appear for higher particle inertia near the location of the maximum of 
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(^'i^p,i) (10 < < 20). As observed for the particle velocity rms, the present model as well 
as the one by Simonin et ai— lead to an overprediction of the streamwise component while 
a better agreement is obtained with the third model [Eq. (l2T|) with Di = 0]. Due to the 
asymmetry of the fluid seen-particle velocity correlations, the two non-diagonal components 
are presented in figures M and M The evolution of {u'^v'^ ^) function of the particle 
inertia is well reproduced. Nonetheless, the stochastic models generally underestimate the 
DNS results. A better agreement is obtained for the other non-diagonal component. As 
noted for the majority of the statistics presented in this study, the model by Simonin et al- 
provides an acceptable but less accurate prediction. 

V. CONCLUSION 

We present in this study a stochastic model for estimating the fiuid velocity experienced 
by small solid particles in a non-homogeneous turbulent flow. In the first part, a new 
stochastic model which is compatible with the limits of the transport equation of the drift 
velocity for low and high particle inertia has been derived. From this compatibility criterion, 
it has also been shown that some previously proposed stochastic models should not be able 
to refiect accurately the inertia eff'ect on the fiuid velocity seen by high inertia particles. 

In the second part of this study, the accuracy of the present stochastic equation has 
been evaluated. Since no models exist to determine the drift and diffusion parameters, 
appearing in the SDE, for non-homogeneous turbulence, they have been deduced from a 
method similar to the one proposed by Pope^. The obtained results show that the drift 
and diffusion matrices are highly space-dependent and anisotropic. Using these values, 
stochastic simulations of a gas-solid channel flow have been conducted and compared to 
DNS. The stochastically predicted data have been also compared to those obtained with the 
model proposed by Simonin et al— and to a simpler one. Using the compatibility criterion 
presented in the first part, the former model should not be able to reproduce the dynamics 
of high particle inertia while the latter should not be able to reproduce it for small particle 
inertia. 

The three models considered are able to predict with a good accuracy the first and second 
order statistical moments of the particle and fiuid seen velocities. Surprisingly, a good 
accordance has been also noticed for the triple particle velocity correlations. This accordance 
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is mainly qualitative. Nevertheless, some components have been seen to be well predicted 
quantitatively. This clearly demonstrates the capability of Langevin-type models to predict 
accurately and efficiently the interactions between inertial particles and turbulence. 
The accuracy of the results obtained with these three different models diverges principally for 
the particle concentration and the drift velocity. As stated before, a good estimation of these 
quantities is primordial to correctly predict the important process of particle deposition. It 
has been seen that the model proposed by Simonin et al- is not able to predict the increase 
near the wall of the concentration of moderate and high particle inertia. On contrary, 
the simpler model predicts this increase even for the smaller particle inertia whereas the 
concentration should be almost uniform. This clearly shows that this model suffers from a 
spurious drift effect. The new proposed model is the only one which succeeds in predicting 
the good evolution of the particle concentration for the range of particle inertia studied. 
This naturally leads us to consider that it is a good candidate to estimate the turbulence 
seen by inertial particles. In the present paper, we test the proposed model using DNS data 
such as for the fluid-particle covariances in order to not introduce supplementary modeling 
uncertainties. These data are generally not known beforehand. Nonetheless, numerical 
strategies related to RANS-Lagrangian methods can help to overcome this difficulty. For 
instance, one such method can be found in Peirano et al.—, which calculates the fluid-particle 
covariances on the basis of the statistics of a large number of particles. Thus, the model is 
more computationally demanding than Simonin's method, but it is hoped that the benefit 
of better results outweighs the extra cost. 

Another possible way to bypass this difficulty would be to directly model the fluid-particle 
covariances. The simplest existing model is based on the theory developed by Tcheu'^ and 
Hinzo^ [see Simonin et al.-]. Attention must be paid to the properties of the carrier fluid 
flow studied since this model was initially developed for isotropic and stationary turbulence 
under restrictive assumptions. A more sophisticated approach was recently proposed by 
Zaichik et al.-^. Although their model was developed for quasi-homogeneous anistropic 
turbulence, quite accurate predictions of the fluid-particle covariances were obtain for a 
gas-solid non-homogeneous flow.— Moreover, the union of these models with the proposed 
stochastic equation is consistent since both depend on the same set of quantities, i.e. the 
decorrelation time scales and second order statistical moment of the fluid velocity seen by 
particles. 
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We would like to emphasize that the model proposed belongs to a particular class of 
stochastic models. Many other models could certainly reproduce more accurately the in- 
teractions between inertial particles and turbulence. Nonetheless, the difficulty is to find a 
model which is a good compromise between complexity and physical accuracy. Besides, the 
challenge lies also in the specification of the model parameters. One could propose a model 
which is theoretically able to reproduce the different physical aspects of gas-solid flows with 
high fidelity, but if its parameters cannot be estimated accurately, the model will probably 
be less satisfactory than a simpler model whose parameters can be found precisely. One part 
of the physics is in the functional form of the model, the other part is in its parameters. 

Finally, it should be noted that Lagrangian stochastic methods can also be considered for 
predicting the motion of gas bubbles in a liquid. Nevertheless, deformation, coalescence and 
break-up can significantly modify bubbles shape, and consequently, alter the interactions 
of each bubbles with turbulence. This makes the use of a stochastic model more complex 
since its parameters, which are not well known, should be correctly modified during the bub- 
ble tracking. In addition, the use of the point-force approximation could become ambiguous 
when bubble coalescence is strong. However, if we restrict ourselves to small non-deformable 
bubbles and neglect break-up and coalescence, a Lagrangian stochastic method can be re- 
tained. As for solid particles, the difficulty will be to properly estimate the parameters of 
the stochastic model. 
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Tp*" = 1 (a), Tp'" = 5 (b), and r+ = 25 (c). DNS: — . Stochastic simulation 
with A = 5 {uW^^k) /dxk (□) ; A = 5 {<u'^) /dxk (A) ; A = (O). ks 

20 Non-diagonal component of the fluid-particle covariance tensor, {u^v'^^t^, for 
Tp*" = 1 (a), Tp*" = 5 (b), and r:^ = 25 (c). DNS: — . Stochastic simulation 
with Di = d (u'iv'p^^) /dxk (□) ; A = 9 {u'iu[) /dxk (A) ; A = (O). [49 
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FIG. 2: Non-diagonal components of the drift matrix, Gij. Tp = 1 ( — ) ] Tp = 5 { ) ; Tp = 25 
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FIG. 5: Particle concentration, Cp, for r+ = 1 (a), r^^ = 5 (b), and = 25 (c). DNS: — . 
Stochastic simulation with Di = dlu'-v'^j/) /dxk (□) ; A = 9 {u'-u'i^) /dx^ (A) ; Di = {<>). 
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FIG. 9: Particle kinetic shear stress, {v'piv'^^^^^ fo'^ "^p — (^)i '^p — ^ (b), and Tp = 25 (c). DNS 
— . Stochastic simulation with Di = d (u[v' ,\ /dx^ (□) ; A = d {u'^u'y.) /dxk (A) ; = (O). 
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FIG. 10: Triple particle velocity correlation, j'Wp^i^p i), for Tjj^ = 1 (a), = 5 (b), and = 25 
(c). DNS: — . Stochastic simulation with Di = dl^u'^v'^f}^ /dx^^ (□) ; Di = d {u[u'^) /dxk (A) ; 
A = (O). 
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FIG. 11: Triple particle velocity correlation, (^'p,2^p,2''^p,2)) = 1 (*)) = 5 (b), and Tp = 25 

(c). DNS: — . Stochastic simulation with Di = d(^u[v'p^ /dxk (□) ; A = 9 (u'-u'f^) /dxk (A) ; 
A = (O). 
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FIG. 12: Triple particle velocity correlation, (t'p^2^p,i''^p,i)' '^p — ^ '^p — ^ (b), and Tp = 25 
(c). DNS: — . Stochastic simulation with Di = d(u[v'pf}j /dxk (□) ; A = 9 (u'^u'f^) /dxk (A) ; 
A = (O). 
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FIG. 13: Triple particle velocity correlation, (^'p^2'^p,i^p 2)1 ^r Tj^ = 1 (a), = 5 (b), and = 25 
(c). DNS: — . Stochastic simulation with A = d (^u'^Vp f^^ /dx^ (□) ; A = d {u[u'^) /dxk (A) ; 
A = (O). 
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FIG. 14: Triple particle velocity correlation, (t'p^2^p,3''^p,3)' '^p — ^ i^)^ '^p = ^ (b), and Tp = 25 
(c). DNS: — . Stochastic simulation with Di = dl^u'^v'^^ /dx^ (□) ; A = d {u[u'^) /dxk (A) ; 
A = (O). 
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FIG. 15: Streamwise drift velocity, (n'l), for = 1 (a), = 5 (b), and Tp = 25 (c). DNS: 
Stochastic simulation with A = d\^Wpk) /^^k (□) ; A = d {u[u'j^) /dx^ (A) ; A = {O). 
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FIG. 16: Wall-normal drift velocity, {u'^) , for t+ = 1 (a), t+ = 5 (b), and t+ = 25 (c). DNS: — . 
Stochastic simulation with Di = dlu[v'pj\ /dxk (□) ; A = d {u[u'j^) /dx^ (A) ; = {O). 
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FIG. 17: Diagonal component of the fluid-particle covariance tensor, (u^v'^^i^^ for r + = 1 (a), 
Tp*" = 5 (b), and r+ = 25 (c). DNS: — . Stochastic simulation with Di = d(u[v'^j^ /dx}^ (□) ; 
D, = d{u[u';, /dxk (A) ; A = (O). 
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FIG. 19: Non-diagonal component of the fluid-particle covariance tensor, {u'lv'p^-^, for = \ (a), 
T+ = 5 (b), and r+ = 25 (c). DNS: — . Stochastic simulation with Di = d (u^v'^ / dx^. (□) ; 
D, = d{u[u'f^) /dxk (A) ; A = (O). 
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FIG. 20: Non-diagonal component of the fluid-particle covariance tensor, {u2Vp^i), for Tp = 1 (a), 
Tp = 5 (b), and = 25 (c). DNS: — . Stochastic simulation with Di 

D, = d (v'yj^) /dxk (A) ; A = (O). 
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